home *** CD-ROM | disk | FTP | other *** search
/ BMUG PD-ROM 1996 Fall / BMUG Fall'96 PD-ROM.iso / Education / Math / MathPad 2.4 / Examples / Complex Roots < prev    next >
Text File  |  1994-06-14  |  4KB  |  121 lines
  1. -- This example gives formulas for quadratic and cubic roots and uses the image command to visualize a complex function
  2.  
  3. ------------------- quadratic roots ---------------
  4. -- Algorithm for real parts of roots is from by W.H. Press, S. Teukolsky et al,"Numerical Recipes".
  5. -- Occasionally 4ac << b, so one of the roots is (erroneously) called 0.
  6. -- This formulation avoids the problem.
  7. -- implemented by David Derbes (derbes@uhuru.uchicago.edu) for MathPad
  8.  
  9. -- Given a quadratic of the form
  10.  
  11. --          a*x^2 + b*x + c = 0
  12.  
  13. -- with real coefficients, find the (possibly complex) roots.
  14. -- Roots are given in the form x + iy.
  15. ~
  16. sgn(x) = 1 when x >= 0,-1 otherwise
  17. D = (b*b - 4*a*c)    -- discriminant
  18. x1 = -(b + sgn(b)*sqrt(D))/(2*a) when D >= 0, -b/(2*a) otherwise
  19. x2 = c/x1 when D >= 0, -b/(2*a) otherwise
  20. y1 = 0 when D >= 0, sqrt(-D)/(2*a) otherwise
  21. y2 = -y1
  22. ~
  23.  
  24. ------------------- cubic roots -------------------
  25. -- Tartaglia's & Cardano's formulae for the roots of a cubic
  26. -- from Universal Encyclopaedia of Mathematics
  27. -- implemented by David Derbes for MathPad, 8 Sept 1993
  28. -- Given a cubic of the form
  29. --
  30. --            a0*x^3 + a1*x^2 + a2*x + a3 = 0
  31. -- 
  32. -- with real coefficients, find the (possibly complex) roots.
  33.  
  34. c1 = a1/a0  -- "normalize", i.e. make leading coefficient = 1
  35. c2 = a2/a0
  36. c3 = a3/a0
  37.  
  38. -- discriminant D; if D > 0, one real root; else three real roots
  39. -- (at most two distinct real if D = 0)
  40.  
  41. a = c2 - (c1*c1)/3.0
  42. b = (((2.0*c1*c1*c1) - (9.0*c1*c2))/27.0) + c3
  43. D = (b*b/4.0) + (a*a*a/27.0)
  44.  
  45. numRealRoots = 3 when D < 0, 1 otherwise
  46.  
  47. dHalf = sqrt(abs(D))
  48.  
  49. sgnp = -1 when ((-b/2.0) + dHalf) < 0, 1 otherwise    
  50. sgnq = -1 when ((-b/2.0) - dHalf) < 0, 1 otherwise
  51.  
  52. p = 0.0 when ((-b/2.0) + dHalf) = 0,
  53.      sgnp*(abs((-b/2.0) + dHalf))^(1.0/3.0) otherwise
  54. q = 0.0 when ((-b/2.0) - dHalf) = 0,
  55.      sgnq*(abs((-b/2.0) - dHalf))^(1.0/3.0) otherwise
  56.  
  57. s = (-b/2.0)/sqrt(-a*a*a/27.0);  theta = acos(s)
  58.  
  59. -- roots of the form x + iy
  60.  
  61. x1 = 2.0*p - (c1/3.0) when numRealRoots = 1 and abs(D) < 1.0e-10,
  62.      (p+q) - (c1/3.0) when numRealRoots = 1 and abs(D) > 1.0e-10,
  63.      2.0*sqrt(-a/3.0)*cos(theta/3.0) - (c1/3.0) otherwise
  64.  
  65. x2 = -p - (c1/3.0) when numRealRoots = 1 and abs(D) < 1.0e-10,
  66.      -(p+q)/2.0 - (c1/3.0) when numRealRoots = 1 and abs(D) > 1.0e-10,
  67.      2.0*sqrt(-a/3.0)*cos((theta/3.0) + 120) - c1/3.0 otherwise 
  68.  
  69. x3 = x2 when numRealRoots = 1,      
  70.      2.0*sqrt(-a/3.0)*cos((theta/3.0) + 240) - c1/3.0 otherwise 
  71.  
  72. y1 = 0.0  -- no matter what, must have at least one real root
  73.  
  74. y2 = (p-q)*sqrt(3.0)/2.0 when numRealRoots = 1 and abs(D) > 1.0e-10,
  75.      0.0 otherwise
  76.  
  77. y3 = -y2 when numRealRoots = 1 and abs(D) > 1.0e-10,
  78.      0.0 otherwise
  79.  
  80. root1 := {x1,y1}:;     root2 := {x2,y2}:;   root3 := {x3,y3}:
  81.  
  82. ------- Enter the coefficients for the cubic here --------------
  83. --            a0*x^3 + a1*x^2 + a2*x + a3 = 0
  84.  
  85. a0 = 35;    a1 = 15;    a2 = -5;    a3 = -20
  86.  
  87. root1:{0.76,0.00}
  88. root2:{-0.59,0.64}
  89. root3:{-0.59,-0.64}
  90.  
  91. ----------- Check the solution------------
  92. include ":incl:complex ops"
  93. z(x) = a0*Ccube(x) + a1*Csqr(x) + a2*x + {a3,0} -- The complex cubic
  94.  
  95. -- confirm that z(x) is zero at the roots
  96. z(root1):{0.00,0.00}
  97. z(root2):{0.00,0.00}
  98. z(root3):{0.00,0.00}
  99.  
  100. -------- check the solution graphically
  101. -- define an array that samples points on the surface of:
  102. --    Z = abs(z(x)) vs. X = real(x) , Y = imaginary(x)
  103. -- This surface should dip to zero at the roots. (The sampled surface will dip near zero but in general is not sampled exactly at a root).
  104.  
  105. surf[ix,iy] =  x:=Cscale(ix,iy), Cabs(z(x)) dim[m,m]
  106.  
  107. -- The index for surf[] runs from 1 to m. Scale the index to get real and imaginary parts ranging from Xmin to Xmax and Ymin to Ymax
  108. scale(i,min,max,nsteps) = (i-.5)*(max-min)/nsteps+min
  109. Cscale(ix,iy) = {scale(ix,Xmin,Xmax,m), scale(iy,Ymin,Ymax,m)}
  110.  
  111. m=12;   -- surface is sampled on an m by m grid
  112. Xmin = -1; Xmax = 1
  113. Ymin = -1; Ymax = 1
  114. Zmin =  0; Zmax = 50
  115. image surf                  -- show image of the surface
  116. plot {root1,root2,root3}    -- show locations of roots
  117.  
  118. plot z({X,0})[1]/Zmax     -- plot the real part of z for real x
  119. -- This should be zero at real roots. On the plotted surface, real roots are located along y=0 so the real cubic plotted in this way should pass though its real roots.
  120. plot 0
  121.